Problema 1

Se tiene un oscilador harmónico en dos dimensiones, $(x,y)$ en coordenadas cartesianas,en el que la fuerza está dada por $\vec{F}=-k\vec{r}$ y se tiene una fricción dada por $\vec{f}_{fric}=\begin{cases} -\gamma\vec{v} &, \mbox{si } |\vec{v}|<1 \\ -\mu\ v^{3/2} \hat{v} &, \mbox{si } |\vec{v}|>1 \end{cases}$. Si $m=1, \gamma=0.1 \mbox{ y } \mu=0.2$ resuelva los siguientes ejercicios:

a)Encuentre las ecuaciones de movimiento en coordenadas cartesianas y polares. No es necesario resolver las ecuaciones en ambos sistemas, elija el que más le guste y con una transformación de coordenadas obtenga los resultados para el otro sistema de coordenado.

b)Utilice el método de Runge-Kutta de 4to orden para encontrar la solución a las ecuaciones de movimiento en el inciso anterior con condiciones iniciales $\vec{r}_0=(0,1)\ $ y $\ \vec{v}_0=(cos\theta,sen\theta)$, donde $\theta=\frac{n\pi}{6}$ con $n=0,1,...,11$ .

c)Genere una función que calcule le energía para distintas condiciones iniciales y detenga la simulación cuanto ésta sea menor al 1% de su valor inicial.

d)Dibuje las gráficas del momento angular para los casos anteriores del inciso b).

e)Tome de manera arbitraria una condición inicial y dibuje el espacio fase $(x,p_x)\ $ y $\ (\rho,p_\rho)$.

Solución

a)Por 2a Ley de Newton se tiene que $m\ddot{\vec{r}}=-k\vec{r}+\vec{f}_{fric}=\begin{cases}-k\vec{r} -\gamma\vec{v} &, \mbox{si } |\vec{v}|<1 \\-k\vec{r} -\mu\ v^{3/2} \hat{v} &, \mbox{si } |\vec{v}|>1 \end{cases}$. Entonces:

En Coordenadas Cartesianas:Tomamos que $\vec{r}(t)=(x(t),y(t))=x(t)\hat{i}+y(t)\hat{j}$,con $\hat{i}=(1,0), \hat{j}=(0,1)$. Con esto, la velocidad y la aceleración son símplemente: $\vec{v}=\dot{\vec{r}}=\dot{x}\ \hat{i}+\dot{y}\ \hat{j} \Rightarrow \vec{a}=\ddot{\vec{r}}=\ddot{x}\ \hat{i}+\ddot{y}\ \hat{j}$. Entonces, las ecuaciones de movimiento estan dadas por:

$$\ddot{x}_i=\begin{cases}-\frac{k}{m}x_i-\frac{\gamma}{m}\dot{x}_i &\text{, si }\ (\dot{x}_1^2+\dot{x}_2^2)^{1/2}<1 \\ -\frac{k}{m}x_i-\frac{\mu}{m}(\dot{x}_1^2+\dot{x}_2^2)^{1/4}\dot{x}_i &\text{, si }\ (\dot{x}_1^2+\dot{x}_2^2)^{1/2}>1 \end{cases}, \text{donde } i=1,2\ \text{ y }\ (x_1,x_2)=(x,y). $$

Para resolver la ecuación hacemos el cambio de variable $\begin{cases}x_1=x(t)\\x_2=y(t) \\x_3=\dot{x}(t)\\x_4=\dot{y}(t) \end{cases} \Rightarrow \begin{cases}\dot{x_1}=x_3\\ \dot{x_2}=x_4 \\ \dot{x_3}=\ddot{x}(t)\\ \dot{x_4}=\ddot{y}(t) \end{cases}$.

Por tanto, el sistema se puede escribir como $\dot{\vec{x}}=\vec{g}(t,\vec{x})$, con $\vec{x}=(x_1,x_2,x_3,x_4)$ y $$\vec{g}=\begin{cases}(x_3,x_4,-\frac{k}{m}x_1-\frac{\gamma}{m}x_3,-\frac{k}{m}x_2-\frac{\gamma}{m}x_4) &\text{, si }\ (x_3^2+x_4^2)^{1/2}<1 \\(x_3,x_4,-\frac{k}{m}x_1-\frac{\mu}{m}(x_3^2+x_4^2)^{1/4}x_3,-\frac{k}{m}x_2-\frac{\mu}{m}(x_3^2+x_4^2)^{1/4}x_4) &\text{, si }\ (x_3^2+x_4^2)^{1/2}>1 \end{cases}$$.

En Coordendas Polares Planas:Tomamos que $\vec{r}(t)=r(t)\hat{r}(t)$,con $\hat{r}=(cos\theta(t),sen\theta(t))$. Con esto, la velocidad y la aceleración son: $\vec{v}=\dot{\vec{r}}=\dot{r}\ \hat{r} + r\dot{\theta}\ \hat{\theta} \Rightarrow \vec{a}=\ddot{\vec{r}}=(\ddot{r}-r\dot{\theta}^2)\ \hat{r}+(2\dot{r}\dot{\theta}+r\ddot{\theta})\ \hat{\theta}$, donde $\hat{\theta}=(-sen\theta(t),cos\theta(t))$. En este caso, tenemos una ecuación diferencial en cada dirección coordenada, pero las ecuaciones son distintas. Después de hacer todos los cálculos resulta que:

$$\begin{array}{rl}\text{para }\hat{r}:\quad & \begin{cases} r\dot{\theta}^2-\frac{k}{m}r -\frac{\gamma}{m}\dot{r} & \text{, si }\ (r^2\dot{\theta}^2+\dot{r}^2)^{1/2}<1 \\ r\dot{\theta}^2-\frac{k}{m}r -\frac{\mu}{m}(r^2\dot{\theta}^2++\dot{r}^2)^{1/4} \dot{r} &\text{, si }\ (r^2\dot{\theta}^2+\dot{r}^2)^{1/2}>1 \end{cases} \\ \text{para }\hat{\theta}:\quad & \begin{cases} -2\frac{\dot{r}\dot{\theta}}{r} -\frac{\gamma}{m}\dot{\theta} &\text{, si }\ (r^2\dot{\theta}^2+\dot{r}^2)^{1/2}<1 \\ -2\frac{\dot{r}\dot{\theta}}{r} -\frac{\mu}{m}(r^2\dot{\theta}^2+\dot{r}^2)^{1/4}\dot{\theta} & \text{, si }\ (r^2\dot{\theta}^2+\dot{r}^2)^{1/2}>1 \end{cases} \end{array} $$

Para resolver la ecuación hacemos el cambio de variable $\begin{cases}r_1=r(t)\\r_2=\theta(t) \\r_3=\dot{r}(t)\\r_4=\dot{\theta}(t) \end{cases} \Rightarrow \begin{cases}\dot{r_1}=r_3\\ \dot{r_2}=r_4 \\ \dot{r_3}=\ddot{r}(t)\\ \dot{r_4}=\ddot{\theta}(t) \end{cases}$.

Por tanto, el sistema se puede escribir como $\dot{\vec{r}}=\vec{g}(t,\vec{r})$, con $\vec{r}=(r_1,r_2,r_3,r_4)$ y $$\vec{g}=\begin{cases}(r_3,r_4,r_1 r_4^2-\frac{k}{m}r_1-\frac{\gamma}{m}r_3,-2\frac{r_3 r_4}{r_1}-\frac{\gamma}{m}r_4) &\text{, si }\ (r_1^2 r_4^2+r_3^2)^{1/2}<1 \\(r_3,r_4,r_1 r_4^2 -\frac{k}{m}r_1-\frac{\mu}{m}(r_1^2 r_4^2+r_3^2)^{1/4}r_3,-2\frac{r_3 r_4}{r_1}-\frac{\mu}{m}(r_1^2 r_4^2+r_3^2)^{1/4}r_4) &\text{, si }\ (r_1^2 r_4^2+r_3^2)^{1/2}>1 \end{cases}$$.

b)Utilizamos el método de Runge-Kutta de 4o Orden. Para eso utilizaremos el código que se implementó en las tareas previas.

In [268]:
import numpy as np
import pylab as pl
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
def arg_rk4(funcion,x0,tlist,args=0,h=0.1):
    """
    Funcion de Integracion por el Metodo de Runge-Kutta a 4o orden. Sus argumentos son: la funcion de la ecuacion diferencial,
    el valor de la funcion al tiempo inicial x0,una lista con el tiempo inicial y el final tlist=[ti,tf] y 
    el paso de integracion como argumento opcional. 
    Puede tener una lista de parametros opcionales args que se pasan a la funcion de la ecuacion diferencial.
    Regresa los tiempos y sus respectivos valores de la funcion como una tupla de arrays t,x. 
    Cabe notar que la funcion, el jacobiano y x0 pueden tomar valores vectoriales, los cuales se representan por arrays.
    Tambien, la funcion debe tener ordenados los parámetros de la forma g(t,x).
    """
    if (tlist[1]-tlist[0])<10.0*h:   #cambia el valor de h para garantizar que se hagan almenos 10 iteraciones
        h=(tlist[1]-tlist[0])*0.1
    tiempos=np.arange(tlist[0],tlist[1],h) #regresa un array
    x = np.zeros((len(tiempos),len(x0)))#regresa otro array del ancho de la dimension de x0 y el largo del array de tiempos)
    x[0,:] = x0 #x[tiempo,coordenada],:==todas las coordenadas
    if args == 0:  #no hay argumentos opcionales dados    
        for i in xrange(0,len(tiempos)-1):
            k1=funcion(tiempos[i],x[i,:])
            k2=funcion(tiempos[i]+0.5*h,x[i,:]+0.5*h*k1)
            k3=funcion(tiempos[i]+0.5*h,x[i,:]+0.5*h*k2)
            k4=funcion(tiempos[i]+h,x[i,:]+h*k3)
            x[i+1,:]=x[i,:]+(1/6.0)*h*(k1+2*k2+2*k3+k4)    
    else:
        for i in xrange(0,len(tiempos)-1):
            k1=funcion(tiempos[i],x[i,:],args)
            k2=funcion(tiempos[i]+0.5*h,x[i,:]+0.5*h*k1,args)
            k3=funcion(tiempos[i]+0.5*h,x[i,:]+0.5*h*k2,args)
            k4=funcion(tiempos[i]+h,x[i,:]+h*k3,args)
            x[i+1,:]=x[i,:]+(1/6.0)*h*(k1+2*k2+2*k3+k4)    
    return tiempos,x 
#regresa una tupla, que la hace inmutable a diferencia de una lista y que, a diferencia de un array, 
#cada entrada puede ser de distinto tipo (heterogéneas)

Ahora programamos las ecuaciones de movimiento en los distintos sistemas coordenados.

In [3]:
def mov_cartesianas(t,x_descartes,k=5.,m=1.,gamma=0.1,mu=0.2): #x_descartes=(x1,x2,x3,x4)=(x,y,vx,vy)
    x=x_descartes #asumiremos que este es un array
    v=np.sqrt(x[2]**2+x[3]**2)
    if v<=1:
        dx=[x[2],x[3],-k*x[0]/m-gamma*x[2]/m,-k*x[1]/m-gamma*x[3]/m]
    else:
        dx=[x[2],x[3],-k*x[0]/m-mu*np.sqrt(v)*x[2]/m,-k*x[1]/m-mu*np.sqrt(v)*x[3]/m]
    return np.array(dx) #regresa un array    
In [4]:
def mov_polares(t,x_polares,k=5.,m=1.,gamma=0.1,mu=0.2): #x_polares=(r1,r2,r3,r4)=(r,theta,vr,vtheta)
    r=x_polares #asumiremos que este es un array
    v=np.sqrt(r[0]**2*r[3]**2+r[2]**2)
    if v<=1:
        dr=[r[2],r[3],r[0]*r[3]**2-k*r[0]/m-gamma*r[2]/m,-2*r[2]*r[3]/r[0]-gamma*r[3]/m]
    else:
        dr=[r[2],r[3],r[0]*r[3]**2-k*r[0]/m-mu*np.sqrt(v)*r[2]/m,-2*r[2]*r[3]/r[0]-mu*np.sqrt(v)*r[3]/m]
    return np.array(dr) #regresa un array    

Ahora, utilizamos el metodo de Runge-Kutta para integrar las ecuaciones de movimiento, utilizando condiciones iniciales $\vec{r}_0=(0,1)\ $ y $\ \vec{v}_0=(cos\theta,sen\theta)$, donde $\theta=\frac{n\pi}{6}$ con $n=0,1,...,11$ .Además de simplemente dibujar todas las soluciones en una o varias gráficas optaremos por animar la solución a fin de que el resultado sea más ilustrativo.

Coordenadas Cartesianas

In [5]:
paso=0.01
tf=10
x0,y0=0,1 #px=(x,y,vx,vy)

fig = plt.figure(figsize=(5,5))
colormap = plt.cm.gist_rainbow
plt.gca().set_color_cycle([colormap(i) for i in np.linspace(0, 0.9, 12)])
sols_cart=[]
tiempos_cart=np.arange(0,tf,paso)
for n in xrange(12):
    theta=np.pi*n/6
    vx0=np.cos(theta)
    vy0=np.sin(theta)
    px=[x0,y0,vx0,vy0] 
    tx,x = arg_rk4(mov_cartesianas, px, [0,tf],h=paso)
    sols_cart.append(x)
    plt.plot(x[:,0],x[:,1],lw="2")
plt.plot(0,1,marker='o',color="black",markersize=10)
plt.title('Solucion en Coordendas Cartesianas')
plt.show()
In [6]:
from JSAnimation import IPython_display
from matplotlib import animation
In [7]:
# Se define el ambiente en el que queremos hacer la animación
fig = plt.figure()
ax = plt.axes(xlim=(-2, 2), ylim=(-2, 2))
lines = [plt.plot([], [],lw=2)[0] for j in range(12)]
points = [plt.plot([], [],marker='o',color="black")[0] for j in range(12)]
T = 24*tf # periodo para generar 2*T cuadros

# Funcion para inicializar cada cuadro de la animacion
def init():    
    for line in lines:
        line.set_data([], [])
    for point in points:
        point.set_data([],[])
    return lines,points

# Esta funcion se llama de manera secuencial para cada elemento k.
step=len(tiempos_cart)/(2*T)
def animate(k):
    for j,line in enumerate(lines):
        xdata=sols_cart[j][:,0]
        ydata=sols_cart[j][:,1]
        line.set_data(xdata[1:(step*k)],ydata[1:(step*k)])#grafica una linea con los puntos previos
    for j,point in enumerate(points):
        xdata=sols_cart[j][:,0]
        ydata=sols_cart[j][:,1]
        point.set_data([xdata[step*k]],ydata[step*k])#grafica sólo el punto actual
    return lines,points

# Se llama a la animacion.  blit=True es para que solo se dibije las partes de la imagen que tienen cambios.
animation.FuncAnimation(fig, animate, init_func=init,frames=2*T, interval=20, blit=True)
Out[7]:


Once Loop Reflect

Coordenadas Polares

En esta caso integraremos en el espacio $(r,\theta)$ y transformaremos la solución de vuelta a coordenadas cartesianas del espacio $(x,y)$. Para esto definimos una función que dado un vector de la forma $(r,\theta,v_r,v_{\theta})$ nos de uno de la forma $(x,y,v_x,v_y)$. Para esto basta recordar que $\vec{v}=\dot{\vec{r}}=\dot{x}\hat{i}+\dot{y}\hat{j}=\dot{r}\hat{r}+r\dot{\theta}\hat{\theta}$

In [8]:
def polares_a_cartesianas(r):
    x=r[:,0]*np.cos(r[:,1])
    y=r[:,0]*np.sin(r[:,1])
    vx=np.cos(r[:,1])*r[:,2]-np.sin(r[:,1])*r[:,0]*r[:,3]
    vy=np.sin(r[:,1])*r[:,2]+np.cos(r[:,1])*r[:,0]*r[:,3]
    x_transf=np.zeros((len(x),4))
    for i in xrange(len(x)):
        x_transf[i,:]=np.array([x[i],y[i],vx[i],vy[i]])
    return x_transf
In [9]:
paso=0.01
tf=10
x0,y0=0,1 #px=(x,y,vx,vy)

fig = plt.figure(figsize=(15,7))
fig.suptitle("Solucion en Coordenadas Polares \n",fontsize=18)
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
colormap = plt.cm.gist_rainbow
plt.gca().set_color_cycle([colormap(i) for i in np.linspace(0, 0.9, 12)])

sols_pol=[]
sols_transf=[]
tiempos_pol=np.arange(0,tf,paso)
for n in xrange(12): #hay problema en 3 y en 9 puesto que vtheta es muy pequeño (no cero)
    theta=np.pi*n/6
    r0=1
    theta0=np.pi/2
    vr0=np.sin(theta)
    vtheta0=-np.cos(theta)/r0
    if n==3 or n==9:
        pr=[r0,theta0,(-1)**((n-3)/2),0]
    else:
        pr=[r0,theta0,vr0,vtheta0] 
    tr,r = arg_rk4(mov_polares, pr, [0,tf],h=paso) #integramos la ec dif
    x_transf=polares_a_cartesianas(r) #transformamos la solucion de vuelta a cartesianas
    
    sols_pol.append(r) #añadimos cada solucion al conjunto total
    sols_transf.append(x_transf)
    ax1.plot(r[:,0],r[:,1],lw="2")
    ax2.plot(x_transf[:,0],x_transf[:,1],lw="2")

ax1.plot(1,np.pi/2,marker='o',color="black",markersize=10)
ax1.set_title('Polares')
ax2.plot(0,1,marker='o',color="black",markersize=10)
ax2.set_title('Transformadas a Cartesianas')
plt.grid()
plt.show()

Podemos ver que la solución es equivalente a la obtenida en coordenadas polares, transforamda a coordenadas cartesianas, al compararla con la gráfica que obtuvimos de la solución originalmente en cartesianas. Para ser más precisos, vamos a comparar los valores de ambas soluciones para ver cuánto difieren.

In [10]:
values=[]
for i in xrange(12):
    for j in xrange(4):
        values.append(max(sols_transf[i][:,j]-sols_cart[i][:,j]))
max(values)
Out[10]:
0.00032676166951273178

Vemos que la diferencia máxima en todas las soluciones es un punto es menor que $0.0005$. Este error lo atribuimos a la presición numérica de los algoritmos. Sin embargo, nos dice que ambas soluciones también son parecidas. Para disminuir este valor se podría optar por un paso de integración más pequeño o un método de integración de mayor orden.

Ahora procedemos a realizar la animación de la solución en el espacio $(r,\theta)$. Cabe recordar que $\theta$ esta medido en radianes. Para esto basta con modificar la función de animación que se tenía antes.

In [11]:
fig = plt.figure()
ax = plt.axes(xlim=(-1, 1.5), ylim=(-30,30 ))
lines = [plt.plot([], [],lw=2)[0] for j in range(12)]
points = [plt.plot([], [],marker='o',color="black")[0] for j in range(12)]
T = 24*tf # periodo para generar 2*T cuadros

# Funcion para inicializar cada cuadro de la animacion
def init():    
    for line in lines:
        line.set_data([], [])
    for point in points:
        point.set_data([],[])
    return lines,points

# Esta funcion se llama de manera secuencial para cada elemento k.
step=len(tiempos_cart)/(2*T)
def animate_pol(k):
    for j,line in enumerate(lines):
        xdata=sols_pol[j][:,0]
        ydata=sols_pol[j][:,1]
        line.set_data(xdata[1:(step*k)],ydata[1:(step*k)])#grafica una linea con los puntos previos
    for j,point in enumerate(points):
        xdata=sols_pol[j][:,0]
        ydata=sols_pol[j][:,1]
        point.set_data([xdata[step*k]],ydata[step*k])#grafica sólo el punto actual
    return lines,points

# Se llama a la animacion.  blit=True es para que solo se dibije las partes de la imagen que tienen cambios.
animation.FuncAnimation(fig, animate_pol, init_func=init,frames=2*T, interval=20, blit=True)
Out[11]:


Once Loop Reflect

c) Ahora programamos una función que nos calcule la energía de la partícula en el sistema a partir de ciertas condiciones iniciales. Por comodidad utilizaremos el método de Runge-Kutta que habíamos programado para realizar la integración de manera iterada y, solo cuando veamos que la energía en el último punto es menor a la 1% de la energia inicial, truncaremos los valores no deseados.

Coordenadas Cartesianas: En este caso, la energía es: $E=\frac{m}{2}(v_x^2+v_y^2)+\frac{k}{2}(x^2+y^2)$

In [12]:
def energia_cart(px,tlist,porcentaje=1,he=0.01,k=5.,m=1.):
    v2_ini=px[2]**2+px[3]**2
    energia_ini=0.5*m*v2_ini+0.5*k*(px[0]**2+px[1]**2)
    energia=np.array([energia_ini])
    tiempos=np.array([tlist[0]])
    dt=tlist[1]-tlist[0]
    while energia[-1]>=0.01*porcentaje*energia_ini:
        tx,x = arg_rk4(mov_cartesianas, px,tlist,h=he) #hacemos una integracion preeliminar
        v2_nueva=x[:,2]**2+x[:,3]**2
        energia_nueva=0.5*m*v2_nueva+0.5*k*(x[:,0]**2+x[:,1]**2) #calculamos el array con todas las nuevas energias
        if energia[-1]>=0.01*porcentaje*energia_ini:
            energia=np.concatenate((energia,np.delete(energia_nueva,0)))
            tiempos=np.concatenate((tiempos,np.delete(tx,0)))
        else:
            for i in xrange(len(energia_nueva)-1):
                if energia_nueva[i+1]>=0.01*energia_ini:
                    energia=np.concatenate((energia,np.array([energia_nueva[i+1]]))) #añadimos las energias nuevas mayores al 1% de la inicial
                    tiempos=np.concatenate((tiempos,np.array([tx[i+1]]))) #añadimos los tiempos correspondientes
                else:
                    break
        px=x[-1,:].tolist()
        tlist=[tiempos[-1],tiempos[-1]+dt] #tomamos los nuevos valores de tiempo y posicion para volver a integrar
    return tiempos,energia
In [13]:
paso=0.01
tf=10
fig = plt.figure(figsize=(15,6.5))
ax = fig.add_subplot(111)
for i in xrange(12):
    t,e=energia_cart(sols_cart[i][0,:].tolist(),[0,tf],porcentaje=1,he=paso)
    ax.plot(t,100*e/e[0])
ax.set_ylabel('E/E0 (%)')
ax.set_xlabel('tiempo(seg)')
plt.grid()
plt.show()

Coordenadas Polares: En este caso, la energía es: $E=\frac{m}{2}(r^2 v_{\theta}^2+v_{r}^2)+\frac{k}{2}r^2=\frac{m}{2}(r^2\dot{\theta}^2+\dot{r}^2)+\frac{k}{2}r^2$

In [14]:
def energia_pol(pr,tlist,porcentaje=1,he=0.01,k=5.,m=1.):
    v2_ini=pr[2]**2+pr[0]**2*pr[3]**2
    energia_ini=0.5*m*v2_ini+0.5*k*pr[0]**2
    energia=np.array([energia_ini])
    tiempos=np.array([tlist[0]])
    dt=tlist[1]-tlist[0]
    while energia[-1]>=0.01*porcentaje*energia_ini:
        tr,r = arg_rk4(mov_polares, pr,tlist,h=he) #hacemos una integracion preeliminar
        v2_nueva=r[:,2]**2+r[:,0]**2*r[:,3]**2
        energia_nueva=0.5*m*v2_nueva+0.5*k*r[:,0]**2 #calculamos el array con todas las nuevas energias
        if energia[-1]>=0.01*porcentaje*energia_ini:
            energia=np.concatenate((energia,np.delete(energia_nueva,0)))
            tiempos=np.concatenate((tiempos,np.delete(tr,0)))
        else:
            for i in xrange(len(energia_nueva)-1):
                if energia_nueva[i+1]>=0.01*energia_ini:
                    energia=np.concatenate((energia,np.array([energia_nueva[i+1]]))) #añadimos las energias nuevas mayores al 1% de la inicial
                    tiempos=np.concatenate((tiempos,np.array([tr[i+1]]))) #añadimos los tiempos correspondientes
                else:
                    break
        pr=r[-1,:].tolist()
        tlist=[tiempos[-1],tiempos[-1]+dt] #tomamos los nuevos valores de tiempo y posicion para volver a integrar
    return tiempos,energia
In [15]:
paso=0.01
tf=10
fig = plt.figure(figsize=(15,6.5))
ax = fig.add_subplot(111)
for i in xrange(12):
    t,e=energia_pol(sols_pol[i][0,:].tolist(),[0,tf],porcentaje=1,he=paso)
    ax.plot(t,100*e/e[0])
ax.set_ylabel('E/E0 (%)')
ax.set_xlabel('tiempo(seg)')
plt.grid()
plt.show()

Observamos que el comportamiento es similar en coordenadas polares que en cartesianas, por lo que podemos concluir que nuestros cálculos son consistentes.

d) En este caso partiremos de las soluciones que calculamos en el inciso b), por lo que no nesecitamos hacer el proceso de integración. Entonces, definiremos una función que tome como entrada el array completo de la solución y regrese el valor del momento angular en torno del eje z para cada valor correspondiente en forma de array.

En coordenadas cartesianas el momento en torno al eje z es: $p_z=m(xv_y-yv_x)$.

En coordenadas polares se tiene que: $p_z=m(rv_\theta)=m(r^2\dot{\theta})$

In [16]:
def momento_cart(sol,m=1.):
    momento=m*(sol[:,0]*sol[:,3]-sol[:,1]*sol[:,2])
    return momento

def momento_pol(sol,m=1.):
    momento=m*sol[:,0]**2*sol[:,3]
    return momento

Ahora realizamos las gráficas del momento angular en base a las soluciones obtenidas en el inciso b).

In [17]:
fig = plt.figure(figsize=(15,15))
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)
for i in xrange(12):
    t_cart=tiempos_cart
    pz_cart=momento_cart(sols_cart[i])
    ax1.plot(t_cart,pz_cart)
    t_pol=tiempos_pol
    pz_pol=momento_pol(sols_pol[i])
    ax2.plot(t_pol,pz_pol)
ax1.set_title('Cartesianas')
ax1.set_ylabel('pz cart')
ax1.set_xlabel('tiempo')

ax2.set_title('Polares')
ax2.set_ylabel('pz pol')
ax2.set_xlabel('tiempo')

plt.grid()
plt.show()

Podemos observar al comparar ambas gráficas que el momento angular es aparentemente el mismo en coordenadas polares que en cartesianas. Esto sólo nos dice que nuestros cálculos son consistentes. Por otro lado, vemos que el momento tiende a cero en todas las soluciones, lo cual se debe simplemente al efecto de la fricción. Además, entre más se aleje la velocidad inicial de una dirección horizontal (donde $p_z$ es máximo) se observa un decrecimiento en el momento angular inicial. En el punto donde la velocidad inicial esta en dirección radial, el momento inicial es cero.

e)Tomaremos de manera arbitraria una condición inicial en coordenadas cartesianas. A partir de esta, al transformarla, obtendremos las condiciones iniciales en coordenadas polares. Con esto graficaremos los espacios fase $(x,p_x)\ $ y $\ (\rho,p_\rho)$. Cabe notar que, dado que $\vec{v}=\dot{r}\hat{r}+r \dot{\theta}\hat{\theta}=\dot{x}\hat{i}+\dot{y}\hat{j}$, se tiene que: $$\begin{cases}\dot{r} & = & \dot{x}cos\theta &+& \dot{y}sin\theta \\r\dot{\theta} &= & \dot{y}cos\theta &-& \dot{x}sin\theta \end{cases}$$

In [18]:
from random import *
In [19]:
tf=25
paso=0.01
m=1.
trayectorias=4

fig = plt.figure(figsize=(15,7))
fig.suptitle("Espacios Fase \n",fontsize=18)
ax1 = fig.add_subplot(1,2,1)
plt.xlabel('$x$', fontsize=25)
plt.ylabel('$p_x$', fontsize=25)
plt.grid()

ax2 = fig.add_subplot(1,2,2)
plt.xlabel('$r$', fontsize=25)
plt.ylabel('$p_r$', fontsize=25)
plt.grid()
colormap = plt.cm.gist_rainbow
plt.gca().set_color_cycle([colormap(i) for i in np.linspace(0, 0.9, trayectorias)])

for i in xrange(trayectorias):
    x0=[uniform(-1.5,1.5),uniform(-1.5,1.5),uniform(-1.5,1.5),uniform(-1.5,1.5)] 
    #tomamos nuestra condicion inicial en coordenadas cartesianas y la transformamos a polares
    r=np.sqrt(x0[0]**2+x0[1]**2)
    theta=np.angle(x0[0]+x0[1]*1j)
    vr=x0[2]*np.cos(theta)+x0[3]*np.sin(theta)
    vtheta=(x0[3]*np.cos(theta)-x0[2]*np.sin(theta))/r
    r0=[r,theta,vr,vtheta]
    #integramos con ambas soluciones iniciales
    tx,x = arg_rk4(mov_cartesianas, x0, [0.,tf],h=paso)
    tr,r = arg_rk4(mov_polares, r0, [0.,tf],h=paso)
    
    #graficamos los espacios fase
    ax1.plot(x[:,0],m*x[:,2],lw="2")
    ax2.plot(r[:,0],m*r[:,2],lw="2")
    #marcamos las condiciones iniciales en el espacio
    ax1.plot(x0[0],m*x0[2],marker='o',color="black",markersize=10)
    ax2.plot(r0[0],m*r0[2],marker='o',color="black",markersize=10)

ax1.set_title('Cartesianas')    
ax2.set_title('Polares')
plt.show()

Problema 2

Suponga que se tiene un sistema como en el problema anterior, pero ahora la partícula tiene carga $q>0$ y en el origen también se encuentra una partícula con carga $q>0$.

a)Sin considerar fricción, ¿el sistema puede tener órbitas cerradas?. De ser así, encuentre algunas de estas órbitas o muestre que el sistema no puede tener este tipo de órbitas. (Apoye sus argumentos con una simulación).

b)Considerando fricción, ¿el sistema puede tener órbitas cerradas?. De ser así, encuentre algunas de estas órbitas o muestre que el sistema no puede tener este tipo de órbitas. (Apoye sus argumentos con una simulación).

Solución

Para este problema solo haremos la integración en coordenadas cartesianas mediante el método de Runge-Kutta de orden 4. Considerando la carga en el origen se tiene que: $$m\ddot{\vec{r}}=-k\vec{r}+\frac{Q}{r^3}\vec{r}+\vec{f}_{fric}=\begin{cases}-k\vec{r}+\frac{Q}{r^3}\vec{r} -\gamma\vec{v} &, \mbox{si } |\vec{v}|<1 \\-k\vec{r}+\frac{Q}{r^3}\vec{r} -\mu\ v^{3/2} \hat{v} &, \mbox{si } |\vec{v}|>1 \end{cases}, \text{ donde }\ Q\equiv \frac{q}{4 \pi \epsilon_0}.$$

Por tanto, el sistema se puede escribir como $\dot{\vec{x}}=\vec{g}(t,\vec{x})$, con $\vec{x}=(x_1,x_2,x_3,x_4)$ y $$\vec{g}=\begin{cases}(x_3,x_4,-\frac{k}{m}x_1+\frac{Q}{(x_1^2+x_2^2)^{3/2}}\frac{x_1}{m}-\frac{\gamma}{m}x_3,-\frac{k}{m}x_2+\frac{Q}{(x_1^2+x_2^2)^{3/2}}\frac{x_2}{m}-\frac{\gamma}{m}x_4) &\text{, si }\ (x_3^2+x_4^2)^{1/2}<1 \\(x_3,x_4,-\frac{k}{m}x_1+\frac{Q}{(x_1^2+x_2^2)^{3/2}}\frac{x_1}{m}-\frac{\mu}{m}(x_3^2+x_4^2)^{1/4}x_3,-\frac{k}{m}x_2+\frac{Q}{(x_1^2+x_2^2)^{3/2}}\frac{x_2}{m}-\frac{\mu}{m}(x_3^2+x_4^2)^{1/4}x_4) &\text{, si }\ (x_3^2+x_4^2)^{1/2}>1 \end{cases}.$$

A continuación programamos las ecuaciones de movimiento, donde la fricción será un parámetro opcional.

In [20]:
def mov_electrico(t,x,args=0,k=5.,m=1.,Q=3.,gamma=0.1,mu=0.2): #x_descartes=(x1,x2,x3,x4)=(x,y,vx,vy)
    #asumiremos que x es un array
    r=np.sqrt(x[0]**2+x[1]**2)
    v=np.sqrt(x[2]**2+x[3]**2)
    if args==0:
        dx=[x[2],x[3],-k*x[0]/m + Q*x[0]/(m*r**3),-k*x[1]/m+Q*x[1]/(m*r**3)]
    else:
        if v<=1:
            dx=[x[2],x[3],-k*x[0]/m + Q*x[0]/(m*r**3)-gamma*x[2]/m,-k*x[1]/m +Q*x[1]/(m*r**3)-gamma*x[3]/m]
        else:
            dx=[x[2],x[3],-k*x[0]/m+Q*x[0]/(m*r**3)-mu*np.sqrt(v)*x[2]/m,-k*x[1]/m+Q*x[1]/(m*r**3)-mu*np.sqrt(v)*x[3]/m]
    return np.array(dx) #regresa un array

a)Considerando que el sistema no tiene fricción, mostraremos que el sistema puede tener orbitas cerradas apoyandonos en la simulación. Para esto nos fijamos primero en la componente radial del movimiento. Deseamos ver donde hay un radio estable en donde exista un balance entre la fuerza atractiva (del oscilador harmónico) y repulsiva (de la fuerza de Coulomb). Entonces, hallaremos el potencial del sistema y buscaremos un punto mínimo de este. Dado que no hay fricción, esperamos que la energía debería mantenerse aproximadamente constante. Vemos que, en este caso:

$$V(r)=V_{Osc.\ Har.}(r)+V_{Coulomb}(r)=\frac{k}{2}r^2 - \frac{Q}{r}$$

Ahora, buscaremos el mínimo del potencial utilizando el método de bisección que programamos en la tarea 2 con la derivada del potencial con respecto a r.

In [219]:
def biseccion(funcion,intervalo,nmax=20,error=1e-5): #(la función, intervalo [x1,x2],no máximo de pasos ,error permitido) 
    """
    Utiliza el método de biparticion para hallar una raíz de una funcion de una variable en el intervalo [x1,x2] 
    dado como una lista de dos valores en el dominio. No importa si x1<x2 o x2<x1. 
    Como argumentos opcionales toma el numero de pasos a realizar, con el fin de evitar un loop infinito, y el error
    máximo permitido. Este error se usa para evaluar el ancho del intervalo en x de la ultima iteracion.
    Esta funcion asume que hay almenos una raiz en el intervalo especificado, por lo que siempre regresa un valor.
    """
    if intervalo[0]<=intervalo[1]: #primero aseguramos que xi<=xd
        xi=intervalo[0] 
        xd=intervalo[1]
    else:
        xi=intervalo[1]
        xd=intervalo[0]     
    n=1
    while n<=nmax:
        xm = 0.5*(xi+xd)
        if funcion(xm)==0 or 0.5*(xd-xi)<error: #si la funcion es cero o los puntos estan muy cerca entre sí
            break
        else:
            n+=1
            if funcion(xi)*funcion(xm)>0: #no hay cambio de signo en la 1a mitad del intervalo
                xi=xm
            else: #si hay cambio de signo en la 1a mitad del intervalo
                xd=xm
    return xm

Ahora definimos la función de la energía potencial junto con su derivada para hallar el mínimo del potencial y graficamos los resultados.

In [22]:
def potencial(r,k=5.,Q=3.):
    return 0.5*k*r**2 + Q/r

def der_potencial(r,k=5.,Q=3.):
    return k*r - Q/r**2
In [235]:
k=5.
Q=3.
m=1.
fig = plt.figure(figsize=(15,5))
ax1 = fig.add_subplot(121)
ax1.set_title('Energia Potencial')
plt.xlabel('r', fontsize=12)
plt.ylabel('V(r)', fontsize=12)
plt.grid()

ax2 = fig.add_subplot(122)
ax2.set_title('Energia Total')
plt.xlabel('r', fontsize=12)
plt.ylabel('E', fontsize=12)
plt.grid()

intervalo=[0.5,10]
minimo=biseccion(der_potencial,intervalo,nmax=100,error=1e-9)
r=np.linspace(0.1,3,100)
ax1.plot(r,0.5*k*r**2,'-.',label="Oscilador")
ax1.plot(r,Q/r,'-.',label="Coulomb")
ax1.plot(r,potencial(r),label="Total")
ax1.plot(minimo,potencial(minimo),marker='o',color="black")
ax1.legend()

t,x=arg_rk4(mov_electrico,[2,0,0,2],[0,50],h=paso)
radio=np.sqrt(x[:,0]**2+x[:,1]**2)
cinetica=0.5*m*(x[:,2]**2+x[:,3]**2)
ax2.plot(r,potencial(r),label="Potencial")
ax2.plot(radio,cinetica+potencial(radio),label="Energia")
ax2.legend()

plt.show()

De las gráficas anteriores, vemos que la energía total de la partícula se mantiene aparentemente constante. Además, esta se mantiene acotada entre 2 radios durante toda la trayectoria.

Ahora, por simplicidad, tomando la condición inicial $(x_0,y_0,v_{x0},v_{y0})=(x_{min},0,0,x_{min})$ buscamos un tiempo en el que la órbita se cierre, regresando al punto inicial. Graficando la trayectoria de la partícula llegamos mediante ensayo y error a un tiempo final estimado de $t_f=77.46$.

In [179]:
tf=77.46
paso=0.001

plt.figure(figsize=(15.,15.))
x0=[minimo,0,0,minimo]
t,x=arg_rk4(mov_electrico,x0,[0,tf],h=paso)
plt.plot(x[:,0],x[:,1])
plt.plot(x[-1,0],x[-1,1],marker='o',color="black")
plt.plot(minimo,0,marker='x',color="black",markersize=20)
plt.xlabel('x(t)')
plt.ylabel('y(t)')
plt.grid()
plt.show()   

Para evaluar qué tan buena es nuestra estimación calculamos la norma de la diferencia del vector posición con respecto a la posición inicial. Deseamos comparar el punto en el que el cuadrado de la norma de estas diferencias sea mínimo. Para evitar comparar con los puntos iniciales de la simulación, comenzamos desde la entrada 1000 del vector de soluciones.

In [218]:
#comenzamos desde la entrada 1000 para que el cuerpo se aleje lo suficiente de la condición inicial.
dif_pos=np.zeros(len(x)-1000)
dif_vel=np.zeros(len(x)-1000)

for i in xrange(len(dif_pos)):
    dif_pos[i]=(x[i+1000,0]-x[0,0])**2+(x[i+1000,1]-x[0,1])**2
for i in xrange(len(dif_vel)):
    dif_vel[i]=(x[i+1000,2]-x[0,2])**2+(x[i+1000,3]-x[0,3])**2    
    
print"len(x)=",len(x),"len(dif_pos)=", len(dif_pos),"len(dif_vel)=", len(dif_vel)
print np.where(dif_pos==min(dif_pos))
print np.where(dif_vel==min(dif_vel))
np.sqrt(dif_pos[76459]),np.sqrt(dif_vel[76459])
len(x)= 77460 len(dif_pos)= 76460 len(dif_vel)= 76460
(array([76459], dtype=int64),)
(array([76459], dtype=int64),)
Out[218]:
(0.00029770798205940437, 0.008906269333306337)

Observamos que la diferencia es mínima en la penúltima entrada de los arrays de diferencias. Además, la norma de estas diferencias es menor a 0.01 en ambos casos. Esto nos dice que la trayectoria sí parece ser cerrada, sobretodo en lo que concierne a la posición.

b)Ahora consideraremos el sistema con la fricción dada en el problema 1. Tomaremos varias condiciones iniciales y veremos el comportamiento de la energía. Dado que las fuerzas solo dependen del radio nuestra posición de inicio estará fija en $(x_0,y_0)=(2,0)$. De manera similar, tomaremos solo velocidades iniciales en la dirección $y$ positiva. Con esto bastará variar las velocidades e integrar cada una de las condiciones iniciales.

In [364]:
tf=50.
fig=plt.figure(figsize=(12.,12.))
plt.xlim([-2.5,2.5]),plt.ylim([-2.5,2.5])
theta=np.linspace(0,2*np.pi,tf/0.1)
xr=minimo*np.cos(theta)
yr=minimo*np.sin(theta)
plt.plot(xr,yr,lw=3,color="black",label="V(r) minimo")
sols_elec=[]
tiempos_elec=[]
for i in xrange(3):
    x0=[2,0,uniform(-5,5),uniform(0,5)]
    t,x=arg_rk4(mov_electrico,x0,[0,tf],args=1,h=0.05)
    sols_elec.append(x)
    tiempos_elec.append(t)
    plt.plot(x[:,0],x[:,1])
    plt.plot(x[-1,0],x[-1,1],color="cyan",marker='o',markersize=10)
plt.plot(sols_elec[0][-1,0],sols_elec[0][-1,1],color="cyan",marker='o',markersize=10,label="Punto Final")
plt.plot(2,0,marker='o',color="black",markersize=5,label="Punto Inicial")
plt.xlabel('x(t)', fontsize=15)
plt.ylabel('y(t)', fontsize=15)
plt.legend()
plt.show()

Como podemos observar, el radio de las trayectorias va disminuyendo progresivamente. Sin embargo, este no tiende a cero. En realidad, las trayectorias tienden al radio en donde el potencial es mínimo, oscilando en torno a este. Debido a la fricción, tanto la componente radial como la tangencial de la velocidad disminuyen progresivamente. Por tanto, en un tiempo infinito, la partícula quedará estacionada en un punto fijo dentro del radio donde el potencial es mínimo.

A continuación graficamos la energía total de las soluciones comparandola con el potencial del sistema.

In [369]:
fig=plt.figure(figsize=(12.,12.))
intervalo=[0.5,10]
r=np.linspace(0.1,3,100)
plt.plot(r,potencial(r),label="Potencial",color="black")
for i in xrange(len(sols_elec)):
    x=sols_elec[i]
    radio=np.sqrt(x[:,0]**2+x[:,1]**2)
    cinetica=0.5*m*(x[:,2]**2+x[:,3]**2)
    plt.plot(radio,cinetica+potencial(radio))
plt.xlabel('r', fontsize=12)
plt.ylabel('E', fontsize=12)
plt.legend()
plt.show()

Podemos observar que la energía total del sistema no es constante en ninguna de las trayectorias y que, efectivamente, decae hasta el mínimo del potencial. Esto es para cualquier trayectoria de las generadas, por lo que podemos asumir que la energía no se conservará siempre que la partícula tenga energía mayor al mínimo del potencial. Por tanto, si una trayectoria empieza fuera del radio mínimo esta siempre decaerá al mismo radio. Incluso, si una partícula inicia su trayectoria en el radio mínimo su velocidad disminuirá progresivamente hasta que esta se quede estacionada. Por lo tanto, podemos concluir que no hay trayectorias cerradas en este sistema, salvo que consideremos partículas con posición constante en el tiempo dentro del radio mínimo.